### Table 1 ###

## Model 1 ##
## Read in the data ##
setwd("")
data2 <- read.csv("womens.csv", header=TRUE)
attach(data2)

## Changing the id2 variable into dummies for FE ##
data2$id2.factor <- as.factor(data2$id2)

## Loading the nlme library to run the GLS model ##
library(nlme)

## Running model 1 from table 1 and summarizing the results ##
model1 <- gls(zFDlaborfemale ~ zFDlogGDPcap_1 + zFDlogGDPcapSQ_1 + zFDage_1 + id2.factor, data = data2, na.action = na.omit, correlation = corAR1(.1, ~year|id2.factor, fixed = FALSE), method = "ML")
summary(model1)

## Model 2 ##
## Running model 2 from table 1 and summarizing the results ##
model2 <- gls(zFDlaborfemale ~ zFDlogGDPcap_1 + zFDlogGDPcapSQ_1 + zFDage_1 + zFDoil_gas_rentPOP_1 + id2, data = data2, na.action = na.omit, correlation = corCAR1(0.2, ~year|id2, fixed = FALSE), method = "ML")
summary(model2)

## Recoding the data to remove the two most influential countries ##
data2cort <- data2[data2$id != "KWT", ]
data2recort <- data2cort[data2cort$id != "SAU", ]

## Changing the id2 variable into dummies for FE ##
data2recort$id2.factor <- as.factor(data2recort$id2)

## Model 3 ##
## Running model 3 from table 1 and summarizing the results ##
model3 <- gls(zFDlaborfemale ~ zFDlogGDPcap_1 + zFDlogGDPcapSQ_1 + zFDage_1 + zFDoil_gas_rentPOP_1 + id2, data = data2recort, na.action = na.omit, correlation = corCAR1(0.2, ~year|id2, fixed = FALSE), method = "ML")
summary(model3)

## Model 4 ##
## Changing the id2 variable into dummies for FE ##
data$id2.factor <- as.factor(data$id2)
data$year.factor <- as.factor(data$year)

## Running model 4 from table 1 and summarizing the results ##
model4 <- lm(zFDlaborfemale ~ zFDlogGDPcap_1 + zFDlogGDPcapSQ_1 + zFDage_1 + zFDoil_gas_rentPOP_1 + id2.factor + year.factor, data = data)
summary(model4)


### Table 2 ###

## Read in the data ##
library(foreign)
data <-read.dta("C:\\Documents and Settings\\anhill\\My Documents\\Gov 2001\\Replication\\Replication data - short.dta")
attach(data)
id


## Loading the lmtest package for coeftest ##
library(lmtest)

## Loading the library package for vcovHC ##
library(sandwich)

## Model 1 ## 
my.lm.2.1 <-lm(zlaborfemaleMOD ~ zlogGDPcap + zlogGDPcap2 + zage15_64 + zme_nafr + zcommie)
my.lm.2.1
summary(my.lm.2.1)

## Finding robust t-values ##
coeftest(my.lm.2.1, vcov=vcovHC(my.lm.2.1, type="HC1"))

## Model 2 ##
my.lm.2.2 <-lm(zlaborfemaleMOD ~ zlogGDPcap + zlogGDPcap2 + zage15_64 + zme_nafr + zcommie + zislam)
my.lm.2.2
summary(my.lm.2.2)

## Finding robust t-values ##
coeftest(my.lm.2.2, vcov=vcovHC(my.lm.2.2, type="HC1"))

## Model 3 ##
my.lm.2.3 <-lm(zlaborfemaleMOD ~ zlogGDPcap + zlogGDPcap2 + zage15_64 + zme_nafr + zcommie + zoil_gas_rentPOP)
my.lm.2.3
summary(my.lm.2.3)

## Finding robust t-values ##
coeftest(my.lm.2.3, vcov=vcovHC(my.lm.2.3, type="HC1"))

## Model 4 ##
my.lm.2.4 <-lm(zlaborfemaleMOD ~ zlogGDPcap + zlogGDPcap2 + zage15_64 + zme_nafr + zcommie + zislam + zoil_gas_rentPOP)
my.lm.2.4
summary(my.lm.2.4)

## Finding robust t-values ##
coeftest(my.lm.2.4, vcov=vcovHC(my.lm.2.4, type="HC1"))


### Table 4 ###

## Read in the data ##
data3 <- read.csv("replic.csv", header=TRUE)
attach(data3)

## Model 1 ##
model5 <- lm(zfemale_seats ~ zlogGDPcap + zme_nafr, data = data3)
summary(model5)

## Finding robust t-values ##
coeftest(model5, vcov=vcovHC(model5, type="HC1"))

## Model 2 ##
model6 <- lm(zfemale_seats ~ zlogGDPcap + zme_nafr + zislam, data = data3)
summary(model6)

## Finding robust t-values ##
coeftest(model6, vcov=vcovHC(model6, type="HC1"))

## Model 3 ##
model7 <- lm(zfemale_seats ~ zlogGDPcap + zme_nafr + zoil_gas_rentPOP, data = data3)
summary(model7)

## Finding robust t-values ##
coeftest(model7, vcov=vcovHC(model7, type="HC1"))

## Model 4 ##
model8 <- lm(zfemale_seats ~ zlogGDPcap + zme_nafr + zislam + zoil_gas_rentPOP, data = data3)
summary(model8)

## Finding robust t-values ##
coeftest(model8, vcov=vcovHC(model8, type="HC1"))

## Model 5 ##
model9 <- lm(zfemale_seats ~ zlogGDPcap + zme_nafr + zislam + zoil_gas_rentPOP + zlaborfemaleMOD, data = data3)
summary(model9)

## Finding robust t-values ##
coeftest(model9, vcov=vcovHC(model9, type="HC1"))

## Model 6 ##
model10 <- lm(zfemale_seats ~ zlogGDPcap + zme_nafr + zislam + zoil_gas_rentPOP + zpolity, data = data3)
summary(model10)

## Finding robust t-values ##
coeftest(model10, vcov=vcovHC(model10, type="HC1"))

## Model 7 ##
model11 <- lm(zfemale_seats ~ zlogGDPcap + zme_nafr + zislam + zoil_gas_rentPOP + zpolity + zPR_MOD, data = data3)
summary(model11)

## Finding robust t-values ##
coeftest(model11, vcov=vcovHC(model11, type="HC1"))

## Model 8 ##
model12 <- lm(zfemale_seats ~ zlogGDPcap + zme_nafr + zislam + zoil_gas_rentPOP + zpolity + zPR_MOD + zclosed_list, data = data3)
summary(model12)

## Finding robust t-values ##
coeftest(model12, vcov=vcovHC(model12, type="HC1"))


### Table 5 ###

## Read in the data ##
data <-read.dta("C:\\Documents and Settings\\anhill\\Desktop\\Replication data - short.dta")
attach(data)

## Model 1 ## 
my.lm.5.1 <-lm(zfemale_min ~ zlogGDPcap + zme_nafr) 
my.lm.5.1
summary(my.lm.5.1)

## Finding robust t-values ##
coeftest(my.lm.5.1, vcov=vcovHC(my.lm.5.1, type="HC1"))

## Model 2 ##
my.lm.5.2 <-lm(zfemale_min ~ zlogGDPcap + zme_nafr + zislam)
my.lm.5.2
summary(my.lm.5.2)

## Finding robust t-values ##
coeftest(my.lm.5.2, vcov=vcovHC(my.lm.5.2, type="HC1"))

## Model 3 ##
my.lm.5.3 <-lm(zfemale_min ~ zlogGDPcap + zme_nafr + zoil_gas_rentPOP)
my.lm.5.3
summary(my.lm.5.3)

## Finding robust t-values ##
coeftest(my.lm.5.3, vcov=vcovHC(my.lm.5.3, type="HC1"))

## Model 4 ##
my.lm.5.4 <-lm(zfemale_min ~ zlogGDPcap + zme_nafr + zislam + zoil_gas_rentPOP)
my.lm.5.4
summary(my.lm.5.4)

## Finding robust t-values ##
coeftest(my.lm.5.4, vcov=vcovHC(my.lm.5.4, type="HC1"))

## Model 5 ##
my.lm.5.5 <-lm(zfemale_min ~ zlogGDPcap + zme_nafr + zislam + zoil_gas_rentPOP + zlaborfemaleMOD)
my.lm.5.5
summary(my.lm.5.5)

## Finding robust t-values ##
coeftest(my.lm.5.5, vcov=vcovHC(my.lm.5.5, type="HC1"))

